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Abstract. Accurate reconstruction of piecewise-smooth functions from a finite number of Fourier coefficients is an important 

C !3 I problem in various applications. This probelm exhibits an inherent inaccuracy, in particular the Gibbs phenomenon, and it 

^*0 I is being intensively investigated during the last decades. Several nonlinear reconstruction methods have been proposed in 

■ the literature, and it is by now well-established that the "classical" convergence order can be completely restored up to the 

^^^' discontinuities. Still, the maximal accuracy of determining the positions of these discontinuities remains an open question. 

C^ ' In this paper we prove that the locations of the jumps (and subsequently the pointwise values of the function) can be 

reconstructed with at least "half the classical accuracy". In particular, we develop a constructive approximation procedure 

which, given the first k Fourier coefficients of a piecewise-C "•" function, recovers the locations of the jumps with accuracy 

~ fe~^ "'■ \ and the values of the function between the jumps with accuracy ~ k~'^ '^^ (similar estimates are obtained 

for the associated jump magnitudes). A key ingredient of the algorithm is to start with the case of a single discontinuity, 

where a modified version of one of the existing algebraic methods (due to K.Eckhoff) may be applied. It turns out that the 

additional orders of smoothness produce highly correlated error terms in the Fourier coefficients, which eventually cancel 

out in the corresponding algebraic equations. To handle more than one jump, we apply a localization procedure via a 

•^^^ I convolution in the Fourier domain, which eventually preserves the accuracy estimates obtained for the single jump. We 

-. ^ ' provide some numerical results which support the theoretical predictions. 

^S^ ■ 1. Introduction 

■ ' Consider the problem of reconstructing a function / : [— 7r,7r] — > M from a finite number of its Fourier coefficients 

def 1 f^ .,., -tkt. 



m 



Ckif) = — / fit) e-^^^ dt, fc = 0, 1, . . . M 

27r J_7r 

It is well-known that for periodic smooth functions, the truncated Fourier series 

M 
def V^ , ,, ^kx 



(N 
> 



OO ■ |J;|=0 

converges to / very fast, subsequently making Fourier analysis very attractive in a vast number of applications. We have by 
the classical Lebesgue lemma (see e.g. [30 ) that 



in 

o, 

O. max |/(x)-5m(/)WI <(3 + lnM)-BM(/) 

— 7r<a:<7r 

where Em (/) is the error of the best uniform approximation to / by trigonometric polynomials of degree at most M. This 
number, in turn, depends on the smoothness of the function. In particular: 

t , • (1) If / is d-times continuously differentiable (including at the endpoints) and /('*)(ir) < R, then (see 1391 Vol.1, Chapter 

C^ I 3, Theorem 13.6]) 

BMU)<Ci-R-M-'^ (1.1) 

(2) If / is analytic, then by classical results of S.Bernstein (see e.g. 30. Chapter IX]) there exist constants C and g < 1 
such that 

EM{f)<C-q'^ (1.2) 

Yet many realistic phenomena exhibit discontinuities, in which case the unknown function / is only piecewise-smooth. As 
a result, the trigonometric polynomial "Sm (/) no longer provides a good approximation to / due to the slow convergence of 
the Fourier series (one of the manifestations of this fact is commonly known as the "Gibbs phenomenon") . It has very serious 
implications, for example when using spectral methods to calculate solutions of PDEs with shocks. Therefore an important 
question arises: "Can such piecewise-smooth functions be reconstructed from their Fourier measurements, with accuracy which 
IS comparable to the 'classical' one (such as l|l.l|l or 1)1. 2|l )"? 
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This problem has received much attention, especially in the last few decades (' [23l[Tn [9l[5l[20l[8l [29l[T4l[T5ll26l l3l l22lll3lll0l l4ll37l 

would be only a partial list). It has long been known that the key problem for Fourier series acceleration is the detection 

of the shock locations. By now it is well-established that classical convergence rates can be restored uniformly up to the 

discontinuities (see e.g. |22p. but the corresponding question for the jump locations themselves is still open. Notice that any 

linear approximation procedure with free (a-priori unknown) jump locations will not be able to achieve accuracy higher than 

-^ - see 1 17 . 
•Jm l - 

Several partial results and conjectures in this direction are known, in particular the following. The concentration method of 

Gelb&Tadmor | 19l [20] recovers the jumps with first order accuracy, and it can be extended to higher orders. Kvernadze | 26 | 

proves that his method can recover jumps of a C^ function with second order accuracy. In [ 171 |7j we have conjectured that 

the locations of the jumps of a piecewise C^ function can be recovered with accuracy k~'^ from its first fc Fourier coefficients 

(a similar conjecture is made in [36 ). Both Eckhoff [14] and Banerjee&Geer [3] made the same conjectures with respect to 

their particular reconstruction methods. We would also like to mention a related but different problem: reconstruction of 

piecewise-smooth functions from point measurements. There, adaptive approximations can achieve asymptotic accuracy fc^"* 

for piecewise C* functions [2| |3T] [27] . 

With this motivation, our main goal in this paper is to arrive at a better understanding of the "optimal", or the "best possible" 
accuracy of reconstruction, especially with respect to the locations and the magnitudes of the jumps. As a means to achieve 
this goal, we develop a reconstruction method which allows for explicit accuracy analysis. Our method is a "hybrid" between a 
Fourier filtering technique which is first applied to localize the jumps, and the algebraic approach of Eckhoff/Kvernadze which 
is used in order to resolve each discontinuity one at a time to a high order of accuracy. It is precisely this "localization" which 
makes the subsequent analysis tractable. 

Our accuracy analysis is "asymptotic" in nature, although we provide the explicit constants at every step. These constants in 
general depend upon various a-priori estimates (such as the minimal distance between the jumps, or the upper bound on the 
jump magnitudes), which are presumably available. See discussion in Section [2] below, in particular 112. 4D . 

Let us now give a brief summary of the main results. 

(1) If a function with a single jump has at least 2(i + 1 continuous derivatives everywhere except the jump, then the 
jump location can be recovered from the first M Fourier coefficients with error at most ~ M~'^~'^ (Theorem 14. 13D . 
In addition, a jump in the i-th derivative can be recovered with error at most ~ M^~'^~^ (Theorem 14.2111 . A key 
observation in the analysis is that the additional orders of smoothness produce highly correlated error terms in the 
Fourier coefficients, which eventually cancel out in the corresponding algebraic equations. 

(2) The localization step does not "destroy" the above accuracy estimates (Theorem [5]2j . Thus, the pointwise values of / 
are recovered with the accuracy ~ M^"*^^ (Theorem 16. l|l up to the jumps. 

(3) Numerical simulations are consistent with the theoretical accuracy predictions (Section [7]). 

By means of this constructive approximation procedure with provable asymptotic convergence properties, we therefore demon- 
strate that the algebraic reconstruction methods for piecewise-smooth data can be at least "half accurate" compared to the 
classical approximation theory for smooth data. 

We provide an overview of the reconstruction procedure in Section |2] For expository reasons, the details of the localization step 
and the analysis of its accuracy are postponed until Section [5] The resolution method of a single jump is presented in Section 
[3] while Section |4] is devoted to proving its asymptotic convergence order. Finally, the accuracy of the whole reconstruction is 
analyzed in Section |6] Some common notations used throughout the paper are summarized below. 

We would like to thank Ch.Fefferman, E.Tadmor and N.Zobin for useful discussions. 



1.1. Notation. 

• N denotes the natural numbers, R - the real numbers, C - the complex numbers. 

• C* denotes the class of smooth functions which are continuously differentiable d times everywhere. C°° is the class of 
smooth functions having continuous derivatives of all orders. 

• Br (z) is the ball of radius r centered at z, and dBr (z) is the boundary of such a ball. 




Figure 2.1. A piecewise-smooth function 



2. The algebraic reconstruction method 



Let us assume that / has K > jump discontinuities {£j}^j. Furthermore, we assume that / £ C* in every segment (fj-i.^j), 
and we denote the associated jump magnitudes at ^j by 



We write the piecewise smooth / as the sum / = * + '!>, where *(ie) is smooth and periodic and ^(ir) is a piecewise polynomial 
of degree d, uniquely determined by {^j} ,{-Ai,j} such that it "absorbs" all the discontinuities of / and its first d derivatives. 
This idea is very old and goes back at least to A.N.Krylov {[25l|4]). Eckhoff derives the following explicit representation for 
*(a;): 



K d 



j = l l-O 



(2.1) 



Vn{x;i,) 



{2^r 



3n+l 



X- ^j 



ij <x< ij + 2-K 



^^' {n + l)\~""- V 27r 

where Vn{x\^j) is understood to be periodically extended to [— 7r,7r] and Bn{x) is the n-th Bernoulli polynomial. For com- 
pleteness, let us dervie the formula for the Fourier coefficients of 'i>(x) (it can also be found in il4j ). 



Lemma 2.1. Let ^[x) be a piecewise polynomial of degree d, with jump discontinuities {£j} -_, and the associated jump 
magnitudes {•A;,j};Zg' "'j • For definiteness, let us assume that co{i) = J^ #(ir)dx = 0. Then 



"^(*) = ^ Zl""''^' ^(^fc)-'-^^^,, 



3^1 



(2.2) 



Proof. One integration by parts yields for k yt 0: 

K 

i=o 



^^\ -27rifc 27rifc i_.^ ^ 

]=0 ^ 



iiE-o.--^^i-(E*^- 



i=o 



3 



and so after d + 1-fold repetition we obtain (recall that # = 0): 



d K 

A key observation is that if ^ is sufficiently smooth, then the contribution of Cfe(^') to Ck{f) is negligible for large k. Therefore, 
for some large enough M one can build from the equations 1)2. 2p an approximate system 

K d 

.,(y)«±^.|^^ fc.M,...,M + <i + l 

j=l 1=0 

Here and in the rest of the paper we use the notation aij — e^^^J . 

In fact, this system (up to a change of variables and the number of equations) lies at the heart of the algebraic reconstruction 
methods of Eckhoff [14], Banerjee&Geer [3] and Kvernadze 26 . Banerjee&Geer solve it for all the parameters at once by least 
squares minimization. Eckhoff and Kvernadze eliminate all the {^4^^^} first, resulting in a system of polynomial equations for 
the {$j}, whose coefficients have nonlinear dependence on the initial data. 

In contrast, we propose to solve this system separately for each ( = (j, beacuse this case reduces to a single polynomial equation 
with respect to f . We achieve this "separation" by filtering the original Fourier coefficients such that only the part related to 
a particular ^j remains. This step requires some a-priori knowledge of the approximate locations of the jumps. Fortunately, 
such an information can easily be obtained by a variety of methods - see Section \5\ 

Let us finish this section by presenting the main steps of the reconstruction. We denote the approximately reconstructed 
parameters with a tilde sign. If not stated otherwise, it is understood that these approximations depend on the index M. It is 
important to note that we distinguish between the actual smoothness of the function / and the reconstruction order. 

Algorithm 2.2. Let f be a piecewise-smooth function with jumps at {^j}"^,, continuously dijferentiable di times between 
the jumps. Fix a reconstruction order to be some nonnegative integer d < di . Let there be given the first M + d + 2 Fourier 
coefficients of f . 

(1) Solve the system l|3.1|l by localization ( Alaonthm \5 . 1\) and resolution (Algorithm I3.g)) . This will provide us with 
approximate values for the parameters \ (j \ and < Ai j | . 

(2) Calculate the sequence 

K d ~ 

cfc(?) = — y^wj=y^^^^ ifci<M 

j = l i=o 
and subsequently recover the approximate Fourier coefficients of the smooth part: 

" def " 

Cfc(*) = Cfc(/)-Cfc(*) |fc|<M 

Take the final approximation to be 

K d 
7=5 + 5= Y^ ck(i)e''"' + J2Y^^l.jM==;i3) (2.3) 

\k\<M j = l 1=0 

The rest of this paper is devoted to providing all the details of the above algorithm and analyzing its accuracy. In particular, 
we will seek estimates of the form 

lU -ij\ < C* {d,di,K) ■ Fi{A,R,G) ■ M"'-'^-'^^^ 

\^i,j-^i.j\ <C**{d,di,K)-F2(A,R,G)- M^'-^''^-'^^^ (2.4) 

\f{x) -f{x)\< C*** {d,d^,K) ■ Fs {A,R,G) ■ MT(''''*i) 

where 

• C*, C**, C*** are some absolute constants depending only on the "size" of the problem; 

• G = G (£i , . . . , (k) represents the geometry of the jump points (such as minimal distance between two adjacent jumps); 

• A = A(|j4o^i| , . . . , |j4di,K|) represents some a-priori bounds on the jump magnitudes, such as lower and upper bounds; 



• iJ is an absolute bound for the Fourier coefficients of the smooth component ^: 

|cfc(*)|<i?.fc-'*-2 (2.5) 

• Fi,F2,F3 and a,P,f are some "simple" functions. 

In the course of our investigation we shall be defining more specific bounds, but it will always be assumed that those can be 
expressed in terms of the above quantities. 

Since we are interested in "asymptotic" estimates, we will in general allow the inequalities ()2.4|l to hold for all M starting from 
some index K* which may be large and depend on the parameters of the problem. However, if a particular bound holds for all 
fc > K* then it will in general hold for k = 1,2,..., K* as well, with some larger multiplicative constants C* . . . , but which 
are harder to compute explicitly. 

3. Resolving a single jump 
Let ui = e~^^. The goal is to recover # from the approximate system of equations 

k ■* 
^^(^) ^ 1;^ E (lijTTT (= ^^(*)) k = M,...,M + d + l (3.1) 

To find 01, we eliminate {j4o, . . . , j4d} from the equations. The result is a single polynomial equation having the exact value 
oj as one of its solutions. In Eckhoff 's paper, this elimination is described in great detail, while here we present only the end 
result. 
Let 

d 

(-) 



iz)'='J2{-iyCy)m,^,z'+^- 



J=0 



Lemma 3.1. The point lo satisfies: 

pf^oj) = Vfc e N 

Proof. The proof is an immediate consequence of Lemma lA. 41 (see Appendix |Aj. ■ 

Since the exact coefficients ruk (and as a result the polynomials p^ (z)) are unknown, we approximate these with the known 
quantities 

rk='2w{tkf+'ckU) 

d+i C3 31 



tiz)'A^J2^-iyCy)r,^,z-+^- 



3=0 
Now we are ready to formulate the procedure of recovering the parameters of a single jump. 

Algorithm 3.2. Let us be given the first M + d + 2 Fourier coefficients of the function f which has a single jump ( G [""'i ■""] • 

(1) Solve the polynomial equation 

gU^) = 

and take ai to be the root which is closest to the unit circle. In Section[^ below, we shall provide the justification for 
this choice. 

(2) The jump magnitudes Aq, . . . , Ad are reconstructed as follows. By p.2|l . the exact values of Aj satisfy 

d 
mfcw"* = y2(ik)''-'Ai Vfc e N (3.4) 

1=0 

5 



We use the approximations r^ ~ mk, cj Ki u and solve the system of linear equations 

d 
r-fcw"* = V^(ifc)'*-' A; k = M,...,M + d (3.5) 



with respect to the unknowns iAif by any one of the standard methods. 



4. Accuracy analysis: a single jump 

Our goal in this section is to analyze Algorithm 13.21 and calculate its accuracy. We shall express all our estimates in terms of 
the index k, keeping in mind that it should be replaced with M to be consistent with the definitions of the previous sections. 

4.1. Accuracy analysis: jump location. We start with the determination of the jump point u. Our strategy will be to 
investigate the polynomial q^{z), and determine the bounds on locations of its roots. We can informally summarize the main 
results as follows: 

(1) Starting from some fc, the roots of g^(z) are "separated" from each other by at least ~ k^^. 

(2) If the function / is continuously differentiable at least di > 2d + 1 times everywhere except at ^, then one of those 
roots deviates from the "true" value oj by at most ~ fc^'*^^. 

We regard gf (z) as a perturbation of P^(z). With this point of view, we shall first describe the roots of p^(z), and then calculate 
the "deviations" due to the difference 

In the subsequent analysis we denote the roots of P^(z) by zj for i = 0, 1, . . . , d, with the convention that Zg = w. Also, 
we denote the roots of g^(z) by k ' . 

It will be convenient to study p'^ (z) in a different coordinate system. For this purpose, consider the following transformation 

of the punctured z-plane: 

u = T{z) = - - 1 z^O 

z 

Then the inverse map is given by 

z = r-i H = -^ « / -1 (4.1) 

ti + 1 

Now we translate the problem into the -u-plane. 
Definition 4.1. For all fc, d G N let 

Claim 4.2. sHu) is a polynomial function. Furthermore, if uq 7^ —1 is a root of sHu), then zo = T~^ ('"o) is a root of pi[u]. 

Therefore it makes sense to study the roots of 5^(u). We denote these roots by cr^ ' ,i = 0, . . . ,d. The observation below is 
an immediate consequence of Lemma l3.ll 

Claim 4.3. s^(0) = 

Therefore we will always take <jg = 0. 

In what follows, we shall break s^(u) into a sum of terms and subsequently apply a perturbation analysis to determine its 
roots. We begin with some simplifications: 



i=0 I i=0 ) 



d+1 d d d+1 



j=o ' i-o i-o j-o 



Now substitute the binomial expansions 

d-l 

fd-h 

m 



Ik + jf-^ = Y^ km.^d-1-m. f \ 

m—o 



^5' 
s=0 



d-l d+1 3 



and get 

stiu) = ^ .-'^, ^ .- ^ (-IF c ; ^) E -' ay-'- r; ') 

/=0 771—0 j—0 s—0 

Now we make a change in indexing according to the following scheme: 

d+l j d+l d+l d d-l d d-m 

j=0 s s=0 j=s 1=0 m=0 771=0 !=0 

and continue: 

d+l d d—m d+l 



^t(«) - E "^ E ^"^ E ('„ 'y-'^' E^-^y (' • y-'- 



S = 771 = 1 = j = s 

Definition 4.4. For all integers t, 5 with < t < d and < 5 < d + 1 let 

j=s 
With this definition, we can write 

d+l d d — m 

4M = E "' E *'" E C^ " )^''-'Al ■ T(d, m + l,s) 

S = 771 = i = 

We will need a technical result. 
Lemma HA. 511 . For all < s < d+ 1 

(1) Ifm + l>s then T{d, m + l,s) ^0 

(2) T{d, 5-1,5) = {-If+Hd + 1 - 5)! (''+^) 

Proof. See Appendix IaI ■ 

The polynomial 5^(u) must therefore be of the form 

d+l 
5^ (u) = 2_^ {°'0,s + ai,sk + • • • + as-i,sk^^^) u^ 

3=1 

where in particular 

as-i,s= {^'^ ^)i''Aoi-l)''+\d+ 1 - s)l.(^^^^) (4.3) 

We break up the polynomial sHu] into a "dominant" and a "perturbation" part: 5^ (u) — 5^ (u) + hf{u] where 

d+l 
sjiu)'^^^ y^^as-i,sk'-'^u' 

JdTi (^•^) 

'^fc(") = / ("0,5 + 11,3*: + \- as-2,sk^~^j u^ 



def ■ 

s=2 



Next we shall see that the dominant component 5^(u) determines the locations of the roots of 5^(u) up to the first order 
accuracy, while the other component h'^{u) is responsible for second-order perturbations of these roots. 

7 



We denote the roots of T^(u) by o-^ ' ,i — 0, . . . ,d with a^ — 0. 
It turns out that 5^(u) can be completely characterized. 

Definition 4.5. For every a > —1 and n = 0,1,2, ... let £„ (x) denote the generalized Laguerre polynomial ([1 Chapter 
22], [34) Chapter 5.2]): 

Lemma 4.6. With the above notations: 
(1) The polynomial 5^(u) satisfies 



n + a\ (— x) 






7fiu) = l7i(ku) (4.5) 



(2) Parthermore, 



Proof. The first part follows from 114. 4D : 

d+l d+l 

k ■ 5^(u) = k ■ y as—i^sk^~ u^ = \ as—i,s {kuY = Sj(fcu) 
s=l s=l 

To prove the second part, we substitute the expression (14.311 into (14. 4D : 

d+l d+l 

,d + U 



5fH = ^a._i,.-u^ = ^(^_ Ji<'>lo(-l)''+l(d + l-5)!( I )u' 

s—l s-1 

d+l 

= -{-ifAo(d + 1)! ^ (^ ^ J _ J ^ = -(-0''^o(<i + l)\4-^\\-u) 



Corollary 4.7. For aH fc £ N, 5^(u*) = i/ and oni?/ i/£^^j-'(-fcu*) = 0. 



|~(M)| 



Lemma 4.8. The numbers \ <^\ ( satisfy the following properties: 

(1) each a^ ' ' is a simple root of s'^ (u); 

(2) Jf'"*' <0 for 1 = 1,2,..., d; 

(S) there exist constants C\ , C2 such that for every fc G N and < i < j < d 



Cifc-^ < 



~(k,d) ~(k,d) 
a- — a 

» 3 



< C2k-^ (4.6) 

Proof. Following Corollary 14.71 we only need to characterize the roots of jC|j~ ^ (— u). By ^34, Chapter 5.2], for every integer 
m > 1 

n\ 
and therefore 



jCj , , (— u) = u- — jC\ ( — u) 



d\ 

•(dTl) 

The polynomials < £„ (e) > form an orthogonal system on the interval (0, 00) (see again T, Chapter 22], 1341 Chapter 

5.2]). Parts ([l]) and ^^ follow immediately. Part ([Sj follows by taking C\ and C2 to be the minimal and the maximal distance 
between the roots of £^ (■"), correspondingly. ■ 

Now we show that ft^(u) perturbes the zeros of 5^(u) by at most ~ fc~^. Since the coefficients ft^ (u) depend linearly on 
j4o . . . ,A.i, we can expect that the bound will depend on the quantity y^,_„ \A.i\. For convenience, let us therefore define 

^*'^=l'max( 1,^1^,1 J 

8 



Lemma 4.9. There exist constants C3,Ki such that for all k > KiA* and for all i = 0, . . . ,d 



~(k,d) (k,d) 
<^i - <^i 



< C3A*k-^ 



(4.7) 



Proof. Our method of proof is based on Rouche's theorem f Theorem IB. 1 [I . We shall define a sequence p{k) = CsA*k ^ (where 
C3 is to be determined) and consider disks of radius p(k) around each one of the roots cr^ . Our goal is to find C3 so that 



ti^o) 



> ?if(ua) for all points tig = (T- ' + p(fc) e^* on the boundaries of these disks. 



In order to bound 



4(^e) 



from below, we shall use Lemma IB. 21 We need to bound from below the first derivative at 

-^ffc d) f'-'ik d^^ 

a^ ' , as well as to bound from above the second derivative in the disk Bj.- 1 I crj ' 
(1) We always have 



du 



ii^) 



~(k.d) d 



-u (r~^(^")) 



d7f(-u;) 



'(k.d) dtU 



~(k,d) 



Now fco-J ' is always a root of sf{u), therefore the value of ^^^(u) 
can write 



~{k,d) 



is independent of k and thus we 



a u 



> Ci — mm 



du 



Since all the roots are simple, this is guaranteed to be a strictly positive bound. 



(2) Now consider a point u* G B^-i I ct^ ' ). Then 
14. 5t and differentiating twice, we get 

d= ~. 



ku* - fc?P'''^ 



< 1 and therefore ku* G Bi I ct^ ] . Using 



d^^^("' 



d-^ ~, 



u=w d^" 



Let 



def 
C5 = max max 



m*gBi 



(~!''") 



d-' -, 
dtii2 



st[w) 



(3) The constants C4 and Cs therefore satisfy the assumptions of Lemma IB. 21 We define Ce = min(l,^j. The 
conclusion is that there exists a constant C7 such that for every function r){k) : N — >■ R satisfying < 7j(fc) < ^ 
we have 



■7;(fc)e'' 



Now we shall bound fe^ ("e) from above. Recall that 

d+l 



> C77j(fc) 



h'l{u)—\ (ao,s + ai.sfc + • • • + as-2,sfc^ ^j u^ 
s=2 

where a^j are some linear functions of Aq, . . . ,Ad. Let i^(fc) : N — > R be any function satisfying < ((^k) < j, and 

consider ug = af'"^"^ + C{k)e^^. By Lemma li^Sl 



~(k,d) 



< C^k ^ and therefore \ue\ < 2 • max(l, C2) k '. But then 

-2 



\hi{ue)\ < |ao,2||«s|^ + (|ao,3| + |ai,3|fc)|us|^ + ---+ (|ao,d+i| + • • • + |ad_i,d+i| fc"*-') \uef+^ < Cg^'fc 
for some constant Cg. 



We set 



C3 



def 2C8 



and let p(fc) — CsA*k ^. We need the inequality p{k) < ^ to be satisfied, and this is obviously possible if 

k>^^A* 
CrCe 

def 



: Ki 



In this case we have shown that 



7i(a^^^^ + p(k), 



> C7p{k) = 2C8A*k- 



and also 



Therefore 



which completes the proof. ■ 

Remark 4.10. We have in fact shown that for each fc > Ki the polynomial si[u) has precisely d+ 1 distinct roots. 

Now we can go back to the original polynomial pfiz) and accurately describe the location of its roots < z^ ' >. Recall from 
Claim 14^21 that z^ — 7 ^ I u^ ' j . Being careful to avoid the singularity a^ — —1 (by choosing large enough fc), we 

now show that the geometry of the roots cr- ' -^ is preserved under 7 ^. In particular, the numbers z- ' remain separated 
from each other (following 114. 6D l, each of them being close (following Lemma 14.911 to one of the numbers 



y(k,d) drf^-l (~(k,d) 



~(k,d) 



(4.8 



The only thing which is different are the constants. 
Lemma 4.11. Let y^ be defined by (14.811 . Then 

(1) there exist constants Cg, Cio, K2 such that for all fc > K2 and < i < j < d 



{k,d) (k,d) 



< Ciofc-^ 



Cgfc-i < y--- y- 

(2) there exist constants C\\,Ks such that for all fc > K3A* 

< CnA* -fc-^ 

(3) there exist constants Ci2,Ci3,Ki such that for all fc > K^A* and < i < j < d 



(k,d) (k,d) 



Ci2k-'^ < 



(k.d) (k,d) 

z) - z,- 



Proof If fc > 2C2 then 



~(k,d) 



< i (see dMl). It follows that ^ < 



Cifc-i < 



< Ciafc-' 

< 1 and so by (14. 8p 



a^''"^ + 1 



vi ' -y) ' 



< iC2k- 



This proves dU with Cg = Ci, Cio = 4C2 and K2 = 2C2 

l-C'.d) I 

< -L 



If in addition fc > -pr^ then 



~(fc,d) (k,d) 



^ — '■ < ^ and therefore 



af"'^' + 1 



(S!,dl (k,d) 

^i - y] 





~fA) _ 


_ ^(Kd) 




zf 


■^' + 1 


^ik,d) 


+ 1 



> I . It follows from (14.711 that 



2C3 



<iC3k-^ fc > max 2C2, ^,^1^4* 



Ci 



■-K3 



and this proves (|2]) with di = 4C3 and K3 as above. 

Let fc > max I K2 , K3 , \ A* . Using ([ij and ([21) , we have one one hand 



def 
= K4 



■i J 



yi'- y) 

(k,d) _ (k,d) 

(k,d) _ (k,d) 



(k,d) (k,d) 

''I - yi 

Cii , ^11 



4fc 2 



(k,d) (k,d) 



(k,d) _ (k,d) 

i ^T 



< -Ciofc-i 



= Ci3 



and on the other hand 



(k,d) (k,d) 



yik.d) _ y{k,d) 



(k,d) (k,d) 



(k,d) (k,d) 



(k,d) 



(k,d) 



> -Cg k-'^ 



- Oi2 



That proves ^. 



Remaining in the z-plane, we now turn to investigate 9^(2) and its roots <k ' >. Recall that we consider gi{z) to be a 
"perturbation" of p^(z) by another polynomial ei{z), i.e. 

qiiz)=pi{z) + eiiz) 

The coefficients of ef (z) depend on the Fourier coefficients of the "smooth part" of our pieceiwise-smooth function /. It turns 
out that in the general setting, the coefficients of e^(z) are large compared to those of p^ (z) and therefore the perturbations 
of the roots are large too. If, however, there is enough structure in those coefficients due to additional orders of smoothness, 
then the perturbation of the roots is small. This is the essense of the key Lemma |4. 121 below. 

Recall that / has in fact di > d continuous derivatives everywhere in [— tt, tt] \ {^}, and denote the additional jump magnitudes 
at ^ by Ad+i , ■ ■ ■ , Ad^ . For every I <. di, let 'i>i denote the piecewise polynomial of degree I with jump point ^ and jump 
magnitudes Ao , . . . ,Ai. Then we write 



/ = *d + (*di - *d) + ** 
where ^* is di-times smooth everywhere in [— 7r,7r]. Thus there exists a constant R* such that 

Let us also denote 



A = max 



(4.9) 
(4.10) 




Lemma 4.12. Let di > 2d + 1. Then there exist constants Ci4,Ci5,K5 such that for all k > K5H 

Cii-H-k-^ i = l,2,...,d 



(k,d) 



(k,d) 



< 



C16 ■ H -k- 



Proof. The idea of the proof is the same as in Lemma [4.91 Namely, we shall seek the constants C14, C15 and K5 such that if 



pi(fc) ^ Ci4 ■ H ■ fc-2 and p2{k) ^ Cn ■ H ■ k' 



(k,d) 



such that p^(z) > e^(z) on the boundary of D^ and 



then for i = 0, 1, . . . , d and k > K5H there exist neighborhoods D, 

(k) 



(k) 



of 



• diamD 

• diamC, 



(fc) 
(fc) 



pi(fc) for ■ 

P2(fc). 



l,2,...,d; 



(_k,d) 



(1) In order to show that p^(z) is large in some neighborhood of z^ ' , let us first show that s^(u) is large in some 
neighborhood of <j^ ' . Recall that 5^(u) — 5^(«) + hf{u). We have shown in the proof of Lemma [4.9l that if 7j(fc) is 
any function satisfying < 7j(fc) < ^, then s^(u) > C7?7(fc) everywhere on dB^^^r^ ia^ ' J . Furthermore, in this 
case ft^(u) < Cs,A*k^'^. We now require that 



CsA*k- 



< -^CTr){k) 



which is true if fc > ^^^ — K-iA* . In that case we have 



~{k,d) 



Jk,d) 



This is almost what we want - we would like to have such a bound on the boundary of a neighborhood of a instead 

"^(k d^ 

of IT. . If i = then these values coincide, and so we're done. Otherwise, recall that for k > KiA* we also have 



~(k,d) 



(k,d) 



< C3A*k ^. So in order to make sure that a ' belongs to B^ri^^ {a ' ), we just require that 

T]{k) > CzA*k-'^. 

We have thus shown the following: 
(a) For every function < r]{k') < ^ and for every fc > K\A* , the following bound holds for every u on the boundary 
of a neighborhood of Cg — of diameter 27j(fc): 



kfcH >-Cir){k') 



(4.11) 



(b) For every fc > K\A* and 7j(fc) as above, which additionally satisfies 7j(fc) > CzA*k^'^, the above bound holds for 
every u on the boundary of a neighborhood of a. ' -^ of the same diameter 2r7(fc), for every i = 0, 1, . . . , d. 
(2) We can now show that similar bounds hold for pf (z)- Again, only the constants will be different. The map 7 ^ l|4.ip . 
being a Mobius transformation, maps B,j(fc) I ct^ J to a circular neighborhood of z^ (which is not necessarily 



centered at z^ ) . Let u* G ^Tfik) \"i\- Now 



u* -a) • ' 



< Cefc-i and also -^ < fff '''' 



< 0. Therefore if 



k > 2(C2 + Co) then iR(u*) > -| and so |u* + 1| > |. On the other hand, in this case |u* + 1| < 2. 

Now let ui and U2 be two points in the u-plane, such that |ui — u^] ~ r and ^ < |-ui| , |-U2| < 2. They are mapped 
to the z-plane such that 

T ftl III 

< 47- 



Ul + 1 U2 + 1 

Recalling (14.1111 and 114.211 , we conclude: 
(a) For every function < ?7(fc) < —fi- and every fc > max (Ki , 2 (C2 + Ce)) j4* , there exists a circular neighborhood of 



w having diameter between -^^^ and 87j(fc), such that the magnitude of p^(z) on the boundary of this neighborhood 
satisfies 

|pf(z)| = Iz-^+i] \sf{u)\ >2-'^-^C7v{k]^Ciev{k) 
(b) For every fc > K'6j4* and 7j(fc) as above, which additionally satisfies 7j(fc) > C3j4*fc~^, the above bound holds for 

(k d) 

every z on the boundary of a neighborhood of z^ of the same diameter as above, for every i = 0,1, . . . ,d. 
(3) Once we have the lower bound for p^(z) on circles of diameter at most Srj (fc) < ^x^ containing z^ , let us now 
establish an upper bound for e^(z) on these circles. Let z* belong to such a circle. On one hand, its distance 



Ak.d) 



,(k.d) 



< ^ for all fc > K^A*. Therefore 



from z^"''"'' is at most — r^ . On the other hand, by Lemma 14.111 

Zt — w| < ^ 1^ ■ Denote C17 — 8Cq + C13 and let zg — (jJ + C{k)e^^ where ^(fc) is some function satisfying 
< C{k) < — r^. Our goal now is to find a uniform upper bound for e^(ze) . 
(a) By 114.911 we have 

Tk-TUk^ 27r(tfc)'*+l {cfc(/) - Cfc(*d)} = 27r(tfc)'*+l {ck{i■d^ ) - Cfe(*d) + Cfc(**)} 



Therefore 



: 27r(ifc) 



d+l 



d+l 



di 

E 



27r ^-^ (ifc)'+l 

l=d+l 



■cj,(**) 



^"■^fl^^'.I^'j^^ 



4M=Y.^-iY{"y) 



Ad+i 



j=o 

dl —d 



di —d 



+ Sk+j 



d+l-j 



def 
= «*. 



1=1 j=o ^ -^-^ i=0 



.^d+i-i 



def 



def 



(b) On one hand, we have the bound H4.10|l . On the other hand, \ze\ < 2 and therefore 

CigR* 



|Afc(ze)| < Ci82''+i • 27r • fc''+^ |cfc (*)| < 



^di —d+l 



for some C19. 



(4.12) 



(c) Now we need to estimate A^ (zg). First 

= 0,'^+^-^ + (d + 1 - j)w'^-^(:{k) e'^ +aj{k) 
where |aj(fc)| < C2oC^ (^) ^°^ some constant C2o- Furthermore, using the estimate of Lemma lA.3l we have 



j=0 



H<C2i-fc-''-'-i 

di — d 






M<C22-fc-<*-' 



di-d d+1 



(fc) 



^ V ' 

\-\<C23-k-H^(k) 

for all k > K? where K7 is an explicit constant (see Lemma FA. 3 II . Therefore 

\Ak {zg)\ < A** X {C24 • fc-''-' + C25C(fc)fc"''"' + C26 • fc"'C' (fc)} 
Combining all the above estimates we therefore have for k > max (K'4, K7] A* 

HM\ < -- X {^ - S^«^) - '-fe (^)) + 5^ (4.13) 

(4) We can finally compare p^ (z) and e^ (z) . Let k > max(K'4, Ky, Ke) j4* and consider two cases. 

def 

(a) Suppose z^ 7^ w. Set p{k) = ^X where C14 is to be determined, and suppose also that 

C3A*k-^ <p{k)<— (4.14) 

k 

We have shown above that there exists a neighborhood D^ containing zj of diameter at most 8/3 (fc) — 
Ci4- H ■ k-^ such that for every z* G dDf''^^ we have |pf(z*)| > Ci6p(fc) = ^^^^ ■ 0° t^e other hand, for 
every such z* we have by H4.13|l 



hf(^*)|<^"x(S^ + SrP« + ^p=W 



Cigi?* 



\^J,d+2 fcd+l'^^ ' k J fcdl-d+1 

^** X (C24 + C2BC6 + C26C|) + C19R* C27 

Therefore we must choose C14 and k for which the condition |l^j* — > ^'^' ,2 ^ is satisfied, together with 

H4.14|l . For example: 

def /8C27 „^ 

Ci4 = max — — , 8C3 

( C'l4\ 

k > max Ki , X H 

\ SCe / 
"^ ^ ' 

def 



In this case, gf (z) has a simple zero k^ in D^ so that 



,(fc,d) (fc,d) 



<Cl4(A*+^"+H*)fc- 



(b) Now consider the case Zg = w. Set jo(fc) = A^, ^ where Cib is to be determined. We again require that 
jo(fc) < —fi-. We have shown that whenever k > K^A* , there exists a neighborhood Dp containing w such that 
for every z* G aD^*"' we have Ipf (z*)| > Ci6p(fc) = '^g|SS^ ' ^"^ ^'^^ °*^^'' '^^"^^ ^y 114^1311 for every such z* 



we have 



A** X (C24 + C25C6 + C26C|) + CigR* C27 {A" + R* 



fcd+a fcd+2 

So we require '^^|g^|^ > ^^''^^+2'^ together with f^i^ < ^ ■ This is possible for example when 

8C27 



Cl5 



C16 



fc > KbH > 



CisH\ d+i 
SCe 



Thus we have completed the proof of Lemma 14.121 



We can finally combine everything and prove the main result of this section. 

Theorem 4.13. Let f have di > 2ci + 1 continuous derivatives everywhere m [— 7r,7r] \ {$}. Let gf(z) be as defined m ()3.3p 



and let 



{"i"'}../' 



enote its roots, such that 



_(k,d) 



< 



Let {<^i}j_i denote the roots of the Laguerre polynomial 

jC^ , such that |(^i| < . . .\4>d\- Let y^ = u and y^ ' ^ = 7 ^ (~k) •^°'" ^ = li • • • i <i (^^^ ()4.8|l ). Then there exist constants 
Cg,Ci5,C2s o.i^d Kg such that for every k > KgH the following statemenets are true: 



(1) The numbers ly^''>lie on the ray Ooj, so that 

Cgk-'^ < 



(k,d) 



(k,d) (k,d) 

vi'- y) 



> 1, and: 

<i < j <d 



(2) Each of the numbers 



AKd) 



{"•'il 



IS close to some y^ 



_(k,d) (k,d) 



{k,d) . 



<C2s-H-k- 



(S) The smallest Kg is very close to w: 

(4) Alaorithm \3.Si vrovides an approximation for cu which is accurate up to order fc-Cii+s) 



Proof. We have already proved Q (for k > K2) and ^ (for k > K3A*) - see Lemma STTl] (0) follows from Lemma l47l2l and 
Lemma r4.11l bv choosing C28 = C'14 + Cii and k > K^H. In order to prove Q, we need to show that no root /c^ is closer 
to the unit circle than /Cg ' . From geometric considerations (see Figure 14.11 on page I15|l , it is sufficient to require that 



,(*=."*) 



(k,d) 



which is true whenever 



Therefore we choose 



<C28-H ■ fc-^ < -Cgfc"^ < - min 
2 2 jT^i 

2C28H 



(k,d) _ (k,d) 



k> 



Cg 



Kg = max [ K2,K3,K5, 



2C2s \ 
Cg J 



4.2. Accuracy analysis: jump magnitudes. Suppose that di > 2d + 1 and let k > KgH so that our algorithm gives an 
approximation a)(fc) with error at most CibH ■ k~'^~^, in accordance with Theorem 14.131 Our goal is to analyze the accuracy 
of calculating the approximate jump magnitudes A, , given by the solution of the linear system (13. 5D . For convenience, we 
denote 



Bi - '' 



^'Ad-l 
i^A^:. 




Figure 4.1. The geometry of j/cf '"^H , ji/i*''''''} , |zi'°''*-'|- The superscripts (k,d) are omitted. The 
picture on the right shrinks towards the unit circle as fc — > oo. 



We consider only the case of exactly d + 1 equations. Thus we can write this system in the following form: 



rfcw 



(k) 



rfc+dw^^ "^ 



Vfc X 



where Vk is the (d + 1) x (d + 1) system matrix 



Vk 



1 k 

1 (fc + 1) 

1 (k + d) 



5(fc) 
^0 



5(fc) 



(fc + 1)'' 
(k + dY 



(4.15) 



By 13. 4D , the "true" coefficients Bj satisfy 



(k) def 



Our goal is to estimate the error s- — Bj — B- ioi j — 0,1,. . . ,d. Let 



rukOJ 

5(fc) 



14 X 



Bo 



Bd 



(k) def k-j -k-j 



Then subtracting H4.16|l from H4.15|l gives 



-0 



Jfc) 



Vo 

(k) 

Vi 



(k) 
"Id 



(4.16) 



(4.17) 



This is the key relation of this section. In order to estimate the magnitude of e ■ , we shall first write out explicit expansions 

(k) —1 

for the quantities tj , and then investigate how these expansions are transformed when multiplied by the matrix Vj. . Our 
analysis will show that the special combination of the structures of both this matrix and the expansion coefficients results in 
remarkable cancellations. 

Let us start by investigating the structure of the matrix Vk- 



Definition 4.14. Let Sj. d denote the (d + 1) x (d + 1) square matrix with entries: 



Example 4.15. For d = 4, we have 



Sk,'t 



/ 1 -fc fc2 -fc3 fc4 \ 

1 -2fc 3fc2 -4fc3 

1 -3fc 6fc2 

1 -4fc 

1 y 



Definition 4.16. For every fc £ N let the symbol Vk denote the following 1 X (d + 1) row vector 



Vk 



[ 1 fc 



With this definition, we have 



14 = 



Vk 

Vk+l 

Vk+d 



(4.18) 



Lemma 4.17. Let fc > 0, then 



V^^ = Sk,d X V-Q-i 



(4.19) 



Proof. Let 1 < m < d + 1 and < t < d. The m-th entry of the vector Vk+t X Sk,d equals to 

771—1 

m — 1 



(vk+t X Sk.d)^ = Ec' + ^Yi-^r-'-' i^_ , _ ;) = t"--' 



and therefore 



SaJqI then follows from JAJEl . 



Vk+t X Sk,d = vt 



(k) 

Now we would like to expand r]- . We can obviously assume the equality 



W(k-) = w ■ 



"(k) 



a{k) 
fcd+2 



Now we estimate wjj.) by Proposition IB. 3l as follows; 



kd+2 J 



where fc is large enough so that 



such that \a{k)\<CiBH 



-k—j 



a(fc)\-(*+J') -fc-7 /, a(fc)a,-i\ "■' _fc_, /, ,, .,a(fc)w-i „,, ., 

-^-41 ^u '' n 1+ \ ,,„ I =t^ '' J l-(fc + j) \;,„ +i;i(fc,j) 



fcd+2 



fcd+2 



is satisfied, and 



oc(k)tj ^ 3 

i.d+2 ^ fc+i+2 

iH. (fc,,)i < (fe+^)(;+^- + i)^-W'^-; < ^^^ . ^.,-.d-s 

2fc2(d+2) f 1 _ ^(fc)c.-'Cfe+i + 2)\ 



Obviously, |rfc| < C30 ■ H ■ k'^. Now by l|4.12|l . we have 
/ di — d 



(k) 



Ad+i 



"—k—j —k—j 



l-l 

dl — d 



-Ik+j 



{^{k + J)y 



mk+, + c.*+^- . J2 , T^\..l + ^^+i I <-"'"' ( 1 - (^+^>W""' ) - m,+,a,-'=-^ + r,+,a.-'=-^/ei (fc) 



fcd+2 



l(fc) 



dj — d 



-4d+! 



^ y B, (fc +,y+^ + V ^^±L^ + H2 (fc 



where |H2(fc,i)| < C31 • H^fc-'*-^ and |^(fc)| < C32 • H. 



Therefore we can write 



(k) 



„W 



(k) 



: P (fc) Bo 



k 
j.d+2 



k+j 
{.d+2 



k+d 
kd+2 



-P(k)Bi 



k' 
■p+2 



(k+j)' 
J.d+2 



(k+d)' 
J.d+2 



-P{k)Bi 



k'^+'^ 
■p+2 



(fc+i)''+ 

J.d + 2 



(k+d)'^^ 
J.d + 2 



-4d+i 



fc+i 



1 

k+d 



Ad+l 






(fc+i)' 






1 
t,d+l 



(fc+j)<* 



(4.20) 



H2(fc,0) 
R2{k,j) 
R2{k,d) 



(k+d)' J L (*:+d)''+l 

In light of l)4.17|l . we would now like to examine the action of the matrix Vg^ on the vectors in the right hand side of l|4.20p . 

Lemma 4.18. Let j = 0,1, . . . ,d. 



(1) If 1 = 1,2,..., d then 



(fc + df 



Vo X 







(2) Otherwise, there exists a function R3 : {0, 1, . . . , d} — > M such that 



(k + j) 



d+l 



{k + d) 



d+l 



Vo X 






H3(0) 



Rs (i) 



H3(d) 



(4.21) 



(4.22) 



Proof. Straightforward application of the binomial theorem. 
Lemma 4.19. For j = 1,2, . . . and i = 1, . . . ,d+ 1 denote 

Then there exists a bounded function R4 : {0, 1, . . . , d} X N — > R such that 



1 

ki 

1 


= Vo X 


r ^' 1 

fc3+l 
^d+1 

L k3 + '' J 


1 


H4(0,i 


(k+iy 


Ri{l,J 


1 


fcd+i+i 


(k+d)3 J 


. Ri[d,3 



(4.23) 



Proof. First recall the well-knowilj power series expansion 



43= 



oo 

1 + xy ^-^^ ' ^ j -1 ' 



( 

n— u 
Now let i = 0,1,..., d. 

(1) On one hand, the [I + l)-st entry in the product on the right-hand side of 114. 23D equals to 



def Y^ 
i=0 



-If ^ .', . e 



(2) On the other hand, by Proposition IB. 31 we have for some bounded function K4 : {0, 1, . . . , d} X N 

{k + iy ki {i + iy fc^ I '^ i-i \^J *''*+' I 



Thus 114. 23D is proved. ■ 

It is now easily seen that the multiplication by Vg" "orders up" the vectors in 114. 20D by decreasing powers of k. Further 
multiplication by Sk,d from the left preserves this structure, as is evident from the following calculation. 

Lemma 4.20. Let Cij be arbitrary constants. Then there exist constants ■yij such that 



kl 
'^2,,- 
W + 1 




V 71,, 
w 

72,, 
W + 1 


'-m - 




Td+1,3 



Sk,d X 



Proof. Let i = 1, . . . , d + 1 and consider the i-th entry of the product, say yi 



d+l 



(4.24) 



- - E (^M)., X -^ = Y.^-ur^-' U L .) X w - ^-'^^ 



1=0 



where 7ij = X)z--i(~''')''''^~'(i-ri-il)'^'+iJ' '^^^^ proves the claim. 
We can now prove the main result of this section. 



Theorem 4.21. Assume that di > 2d + 1 and k > KgH, so that by Theorem \4..13\ we have Wck) — w < C15 ■ H ■ k "* ^. 
Then there exist constants C33, Xio sucft that for every k > KioH and i = 0, 1, . . . d t/ie error in determining At is 



A^-At 



< C33 ■ H^ ■ fc'""*"' 



Proof Combine (gjll), ||42QJ, |g2IJ, ^SSi, llOSJI and (g21J- ■ 

5. Localizing the discontinuities 

As we have seen, both the location and the magnitudes of the jump can be reconstructed with high accuracy. The remaining 
ingredient in our method is to divide the initial function into regions containing a single jump, and subsequently apply the 
reconstruction algorithm in each region. 



It can be proven by induction on j, using the identity N [''"'" ] — ('" ^T" )■ 



Our approach is to multiply the initial function / by a "bump" gj which vanishes outside some neighborhood of the j'-th jump. 
This step requires a-priori estimates of the jump positions, which can fortunately be obtained by a variety of methods, for 
example: 

(1) The concentration method of Gelb&Tadmor 1201 : 

(2) The method of partial sums due to Banerjee&Geer [3]; 

(3) Eckhoff's method with order zero (the Prony method). 

All the above methods provide accurate estimates of {fj} up to first order. For definiteness, we present the description of the 
last method and a rigorous proof of its convergence in Appendix [Dl 

Now, the multiplication is implemented as Fourier domain convolution. Because of the Fourier uncertainty principle, the 
Fourier series of our bump will have inifinite support and therefore every practically computable convolution will always be an 
approximation to the exact one. Nevertheless, an error of order at most fc^'^i^^ in the Fourier coefRcients will be "absorbed" 
in the constant jR* H4.10|l and therefore we will still have accurate estimates for the reconstruction of each separate jump. This 
will require us to use bump functions which are C°°. An explicit construction of such a function is provided in Appendix [Cl 
We assume that the following quantities are known a-priori: 

• the lower and upper bounds for the jump magnitudes of order zero: Ji < |j4o,j| < J2', 

• the minimal distance between any two jumps |^i — $j| > J3 > 0; 

• a constant T for which 



27r (tfc) Cfc (/) - ^ >lo. 



fc 



Our localization algorithm can be summarized as follows. 



<T-fc-^ (5.1) 



Algorithm 5.1. Let f be a piecewise-smooth function of order di > 2d + 1 with K jumps {$j}._, and jump magnitudes 
{•^' ,i }; =0 ' d • ^^^ there be given the Fourier coefficients {c^ (/)},, .jt, , where M is large enough (see below). 

(1) Using the a-prion bounds Ji, J2, J3,T, obtain approximate locations of the jumps \ij\ via Alaorithm \D.3[ In partic- 
ular, the error \^j — ^j should not exceed -^ , and this will be possible if M is not smaller than required by Theorem 



3 

El 



(2) For each $j : 

(a) Construct the bump gj centered at ^j with parameters t = 2- -^ and E = J3, according to Appendix\^ Calculate 
its Fourier coefficients in the range k = —3M . . . 3M according to (|C.2[I . 

(b) Now let hj = f ■ gj . For each fc = 0,l,...,M-|-d + l calculate 



2M 



^k^'^\h,)= J2 ^^if)^k-^i9:) (5.2) 

i=-2M 

(c) Use the above approximate Fourier coefficients ck Cij) <^s the input to Alaonthm \3.S\ for reconstructing all the 
parameters of a single jump. 

Theorem 5.2. Alaorithm \5.1\ will produce estimates of all these parameters with the accuracy as stated m Theorems \4. 1 3\ and 
\4.21\ R* being replaced with some other constant R* — R* [R* ,T, J2, Js]- 

Proof. It is clear that the exact function hj — f ■ gj has exactly one jump at £ and jump magnitudes Aqj, . . . , A^^j. In order to 



prove that Theorems 14. 131 and 14.211 can be applied, it is sufficient to show that the error 
j.-(di+2) gy ^.j^g Fourier convolution theorem the exact Fourier coefRcients of hj are equal to 



Ck (hj) - Cfc (hj 



is of the order 



Ck(.hj) = ^ Ci{f)ck-t{gj) 
i— — 00 

while our algorithm approximates these by the truncated convolution ()5.2|l . Let us estimate the convolution tail 

-2M 00 



Ac['^\h,)^ Y^ c,U)ck-i(g:) + Y^ c,Uyk-.{gj) 



On one hand, the Fourier coefficients of / can be bounded using 115.111 : 

\ck{f)\ <C34(72+T)fc-1 

On the other hand, taking a = di + 1 we have by Theorem lC.il 



kfctoOl < 



C35 1 

jdi+l ' Udi+2 



Finally 






C36 {J2 + T) 



i ( J2 + T) sp _J^ 

^•i .... 



+3 - 



< ^!il;^^c(<ii +3,2M) < £^id^i±n..-^.-^ 



jdi+i 



jdi+i 



-M- 



3 i=2M 3 3 

where C(^i9) i^ the Hurwitz zeta function. Therefore Algorithm 13. 21 will produce estimates of ■j^j > and ■{>li,j } with accuracy 
as guaranteed by Theorems 14. 131 and 14.211 where 



rdl+l 



6. Final accuracy 



In this section we are going to calculate the overall accuracy of approximation. Let us therefore suppose that di > 2d+ 1, and 
so using the Fourier coefficients c_2m(/), • • • 1 C2m(/) we have reconstructed the singular part ^(x) with accuracy specified by 
Theorem 15. 21 Recall from 1)2. 3|l that our final approximation is defined by 

7= Y^ (cfc(/) - Cfc(5)) e'*=- +? 

\k\<M 

where 4 — V^ ._ 5"!;— n ■^l,j^l{^> (j)- Intuitively, the approximation error function f — f will look as depicted in Figure [6711 on 
page l2ll - very small almost everywhere except in some shrinking neighborhoods of the jump points. Let y G [— "", tt] \ {fj} ■_, • 
If we take M large enough so that the error estimate of Theorem l4.13l will be less than the distance to the nearest jump \y — ij\, 
then y will lie in the "flat" region of Figure 16.11 on page [21] and the error / (y) — f (y) will be small. This is precisely the 
content of our final theorem. 



Let us denote the "jump-free" region by 



Dr ='[-7r,7r]\ 



IJ Br Hi) 



\J = 1 



Theorem 6.1. Let f : [— 7r,7r] — >■ M have K jump discontinuities {ij}^-^, and let it be di-times continuously differentiable 
between the jumps. Let r > 0. Then for every integer d satisfying 2d+l < d\, there exist explicit functions F = F (A, R*\ , G = 
G {K,R*,r\ depending on all the a-priori bounds such that for all M y G Alaonthm \2.S\ reconstructs the locations and the 
magnitudes of the jumps with accuracy provided by Theorem \5.Si and with the pomtwise accuracy 



\f{y)-f{y)\ <F-M- 



3/ eDr 



Proof. Define 

|fc|<M 

We write the overall approximation error as 

\ny) - f{y)\ < \f{y) - fM{y)\ + Ifmiy) - f{y)\ (e.i) 

Let us examine the two terms on the right-hand side separately. 

According to our previous notation, * = / — '!> is a d-times continuously differentiable everywhere function. The term \fM — f\ 
is easily seen to be the usual Fourier truncation error of *, since 



\fM{y) - f{y)\ ^ 



J2 (cfe(/)-Cfc(*))e^*=«-^*(t,)-/(t/) 



|fc|<M 



)J CkWe"'y-^{y) 


— 


)J Cfe(*)e'*=« 


k\<M 




\k\>M 




Figure 6.1. The approximation error 



Now recall that the Fourier coefficients of ^ G C* are bounded by 112.511 , therefore 

I/m(2/)-/(j/)| < Css ■ R ■ M-"-"- 



(6.2) 



Let = # — # denote the "singular error function" (see Figure l67l1 on page[2TJ. Then the second term can be written as: 



Ifiy) - fM(y)\ 



Y^ (cfc(*)-Cfc(¥))e'*=''+(?(y)-#(l/)) 



\k\<M 



Y^ Cfe(e)e^*^-e(y) 



|fc|<M 



Write 



5 =£,+a(M) |a(M)| < Fc [a* , A" , R*) ■ M'^-^ 

Ai^j^Ai^j+Pi{M) Wi{M)\ <Ffi[A*,A**,^) • m'-"*-' 

where Fa and Fp are provided by Theorem l5.2l For every e < r, we define 

def 

Using the formula 112.111 we therefore have 

K d K d 

®(^) ^J2Y^{^'''^' iy'^o)-Ai.jVi(y;i,)} = ^^{(^,,, +A(M)) {v, (yUj) + Ui^a(M)(y)) - A,,,V, (yU,)} 



j = l 1=0 
K d 



j = l 1=0 



^^A(M)Vi(3/;$,) + ^^l,,,f7,,„(M)(3/) 



j = l 1=0 



j = l 1=0 



■-z{v) 



■-W{y) 



and so 



\f{y)-fM{y)\ < 



J2 '^ki^)^'-'"' 

\k\>M 



J2 Ck{W)e'''y-W{y) 

\k\<M 



(6.3) 



The functions Vj belong to C' , and therefore by the well-known estimate (see also HI. HI ), there exist constants 5; such that 



J2 ^kivoe'^y 

\k\>M 



<Si-M- 



and therefore 



Y '^kiZ)^'-''^ 



< C39 -Fp-M- 



(6.4) 



\h\yM 



Let us now investigate W {y). Let L denote an upper bound for the magnitudes of the jumps: 

\Ai^j\<L j = \,...,K; l = 0,l,...,di 
Clearly the functions Ui^^ [y] satisfy: 

(1) Since y G Dr, then \Ui eiy)\ < Cio^ for some absolute constant do- This bound can be obtained by just Taylor- 
expanding the functions Vi (y; ^j + e) at e = 0. In particular for e = a (M) we have 

K d 

\W{y)\<J2J2\-^'-'\ \^lMM)iy)\ < C41 • L • Fc • M-''-^ (6.5) 

j=i 1=0 

(2) In the "no man's land" of length a (M) between (j and ^j, Ui^^ is bounded by C42 • L. Furthermore, as we have just 
seen, in the flat regions Ui^ci(M) is bounded by doFa.M^'^^'^ . Therefore the Fourier coefRcients of W are certainly 
bounded by 

\ck{W)\ <Ci3-L-Fc ■M-'^-^ 
and so 



Y^ ck{w)^'^y 



(6.6) 



|fc|<M 
Combining Jell , (lOll , (lOll , (1^41 , (lasl land (lOl l completes the proof. ■ 

7. Numerical results 

In this section we present results of various numerical simulations whose primary goal is to validate the asymptotic accuracy 
predictions for large M. We have used a straightforward implementation and made no attempt to optimize it further. In 
particular, the Fourier coefficients are assumed to be known with arbitrary precision (this is important for the localization, see 
below) . 

7.1. Recovery of a single jump. Given d, di and M, a piecewise function with one discontinuity is generated according to 
the formula 

dl M 

f[x)^Y^AiVi{x;i)+ Y^ /fce^*^ 
1=0 k=-M 

where the numbers ^ G [— 7r,7r],{j4i G IK};lg and {/^ G '^}jf=-M ^^^ chosen at random, such that /j. ~ fc^'*!^^ ^^^ j_^ _ j^ 
The Fourier coefficients are calculated with the exact formula 

1=0 ^ ' 

These coefficients are then passed to the reconstruction routine for a single jump, of order d. This routine implements Algorithm 
I3.2l in a standard MATLAB environment with double-precision calculations. 

The following experiments were carried out: 

(1) Keeping d and d\ fixed, compare the accuracy of recovering the jump location and all the jump magnitudes for 
different values of M. The results can be seen in Figure [7711 on page [23] We also plot the distribution of roots of the 
corresponding polynomials gfj (z) - compare with Figure [4711 on page llBI 

(2) Keeping di fixed, compare the accuracy for different reconstruction orders d= 1, . . . , di. The results are presented in 
Figure [7721 on page [23] 

(3) Keeping the reconstruction order d fixed, compare the accuracy for different smoothness values d\. The results are 
presented in Figure [7751 on page [23] 

The optimality of d = -j- — 1, as well as the asymptotic order of convergence, are clearly seen to fit the theoretical predictions. 
The instability and eventual breakup of the measured accuracy for large values of M is due to the finite-precision calculations. 



Reconstruction accuracy 
d=3,d,=11 
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(a) Accuracy of reconstruction with d = 3 
and di — 11, as a function of M. 



Roots ot qf^(z), d -11 



(b) The roots of 9j^(z) 

Figure 7.1. Reconstruction of a single jump 
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Figure 7.2. Dependence of the accuracy on the order with fixed smoothness, with increasing M. 
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(B) d = 4 



Figure 7.3. Dependence of the accuracy on the smoothness with fixed order, with increasing M. 
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Figure 7.4. Localization: accuracy of recovering the jump. The predicted accuracy M '^ ^ is drawn for 
comparison. 



7.2. Localization. We have restricted ourselves to the following simplified setting: the function has two jumps at f i = and 
^2 ~ 3, and we localize the jump at the origin by a bump having width ^ around the initial approximation ^i ~ -zn- T^^^ 
explicit formulas for the Fourier coefficients of the bump are derived in Appendix [C] We have used Mathematica in order to 
carry out the computations with arbitrary precision. 

The results can be seen in Figure 17.41 on page [24] Localization convergence can clearly be seen here, although it starts from 
very large coefficients. 

8. Discussion 

In this paper we have demonstrated that nonlinear Fourier reconstruction of piecewise-smooth functions can achieve accuracy 

with asymptotic order of at least half the order of smoothness. As indicated by our theoretical results as well as the numerical 

simulations, a reconstruction method whose order is more than half the order of smoothness becomes less accurate. So it appears 

that the algebraic approach has certain limitations, and the interesting question is whether these limitations are inherent or 

superficial. We hope that our results may provide a clue towards obtaining sharp upper bounds. 

In addition, it seems that Eckhoff's conjecture is false as stated in [14 , namely that the jumps of a piecewise-smooth C* 

function can be reconstructed with accuracy fc^"*^^. Using a method of highest possible order doesn't take into account the 

stiffness of the problem. In fact, it can be shown that the Lipschitz constant of the solution map {cfc(/)}j,_^ — > {^ji ^4;^} 

of order d is proportional to M**. We plan to present these results elsewhere. 

Hopefully, our analysis can be related to the algebraic reconstruction schemes of Kvernadze and Banerjee&Geer as well. 

We would like to point out the connection of the algebraic system 112. 2D as well as the well-known Prony system of equations 

IjD.ip (which plays a central role in many branches of mathematics - see [33] and | 28 | ) to other recent nonlinear reconstruction 

methods in Signal Processing, in particular: finite rate of innovation techniques II35II12I , reconstruction of shapes from moments 

1241 121 and piecewise U-finite moment inversion [6][7]. We therefore hope that our results can be extended to these subjects 

as well. 
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Appendix A. Discrete difference calculus and related results 

In this appendix we provide proofs of several combinatorial auxiliary results. 

Let E denote the discrete "shift" operator in k, i.e. for every function g{k) : R — >■ R we have 

Furthermore, let A denote the discrete difference operator, i.e. A = E — I where I is the identity operator. Then by the 

binomial theorem we have 

d 

A^g(k) = (B-I)S(fc) = (-irJ2^-iy(^)g(k+j) (A.l) 

Lemma A.l. Let p(k) = aofc" + aifc"^^ + • • • + a„ 6e a polynomial of degree n. Then 

A'^p{k) = aon! 
A"+lp(fc) = 

Proof. See e.g. [16]. ■ 



[12: 

[13 

[14; 

[is; 
[is: 

[17] 

[is: 

[19 
[20' 
[21 
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[23: 

[24 

[25: 
[26: 

[2?: 
[2s: 

[29 

[30' 
[31 

[32: 
[33: 

[34: 
[35: 

[36: 

[37] 

[3s: 
[39 



Now assume that g (fc) : R"'" — > K is a given function. Let us perform a change of variable y — t, and define G(y) — g ( — ) ~ 
g{k). With this notation, we have 

Ag(k) = gik + 1) - g{k) = g (^- + 1^ - 9 [-) = (^ [j^] - (^ iv) 
We subsequently define the "dual" operator V as: 

V{G{y)} = G(^j^'^-G{y) 

The operator D has an interesting property of "killing" the lowest-order Taylor coefficient at 0. 

Proposition A. 2. Let H(y) be analytic at y = 0, such that H{y) = hmy^ + hm+iy"^~^^ + . . . . Then for n G N, the function 
D" {^H (j/)} IS analytic at with Taylor expansion 

Proof. The proof is by induction on n. The basis n = is given. Assuming that U {y) = D""^ {H(y)'} is analytic at and 
U[y) = «;:;,+„_ 1 3/'"+"- 1 +«;:;,+„_ 1^"^+"- ' + ..., consider the function ©{[/(^y)} = 0"{//(y)}. Let z = ^, and so for 
\y\ < 1 we have 

z{y) = y[l-y + y^ + ...)=y-y^ + ... 
Making the analytic change of coordinates y ^t z {y) we conclude that the Taylor expansion of U (z (y)) at the origin is 

u (z (t/)) = «;,+„_iz'"+"-i + . . . = «:;,+„_i3/-+"-i + . . . 

That is, the leading coefficient is the same as in the Taylor expansion of U [y). Therefore 

V{U{y)} = U{z)-U{y) 



,,m.+n+l 
■'m+n+lH 



This expansion holds in some neighborhood of the origin. 

Lemma A. 3. Let i,d G N. Then there exist positive constants C45,K'ii such that for all k > Kn 

d 



Proof. Let g(k) ~ jr- It is easy to check (see (jA.ip 'l that 



C45 

fcd+! 



3=0 

The proof of the claim is in two steps. First, we shall develop the expression Ai^d(k) into power series in -^ converging for 

sufficiently large k. Then, based on this representation we shall establish the required bound. 

Let y — ^. According to our notation, let G (y) — g (-) — g {k) and so Ai^d{k) — V''' {G (y)} — F{y). Furthermore, for all 

ji = 0, 1, . . . we have 

1 1 + jy 
k+j=-+j= 

y y 

and so 






-jy 

Substituting G (y) — t/', we conclude that F{y] is a real analytic function of y in the disk |t/| < 4, and so it can be written as 

a converging power series 

00 

p{y) = J2 f'y' 

i-O 

Applying Proposition IA.2I to G (t/) we conclude that fo ~ ■ ■ ■ ~ fi+d-i ~ 0- Therefore the expansion 

no 



^'■'*(*')=ES (^-2) 



k 

i—l+d 



26 



holds for k y d. 

Let us now estimate the magnitude of the coefficients fi . Since HA. 211 is valid for fc = d + 1 , then there exists a constant de 
such that \fi (d + 1)~^| < Cue for alH G N and therefore 

\fr\ <Ci6{d+iy 
But then for arbitrary fc > d + 2 we have 



l^!,d(fc)| 



E 



fl+d+i 



f^l+d+i 



< 



C46 (d + 1)'+'* Y^ (d + 1) 



k' 



i=0 



C46 jd + 1]'+'^ 1 £46_ 

d+2 



def 



Lemma A. 4. Let a/ G C, ao, . . . , an G C and n, fc £ N. Denote bk = w*' • (ao + aifc + . . . ank"'). Then 



^ (-IF (-;>.,,.-- 



j=0 



Proof. Denote 

Then we rewrite the given expression as 
n+l 



p(k) — ao + aifc + . . . ank^' 



^ i-iy C ; >.+,-"+l-^ = ^(-1)^- f + >(fc +,ja. 
j=o j—0 

= (-l)'*+la;"-+*+lA"+lp(fc) 

(Lemmc ^A.ll ^ 

The claim is therefore proved. 



n+fc+l 



Lemma A. 5. Let 



d+l 



^K*.^)=T(-^FOr, ) 






Tfte following statements are true for all s = 0,1,. . . ,d+ 1: 

(1) Ift>s thenT{d,t,s) = 

(S) ^(d,5-l,5) = (-l)<'+l(d + l-5)!(''+i) 



def def 

Proof. Let a = d — s and j6 = d — t. Now 

d+l-3 



,j + s^,d+ Ix 



Hdxs)= Y.i-'y^^r:)Qy^^)'-' 

(d+l)! 



j=0 



Z^^ -^ slj! ri-l-.s1!rd-l- 1 - 1 -.s1! 

i=o 



(i + 5)!(d+ 1 - j -5) 



(i + s) 



^(_ix.-c.y-, 1., (1+i)! ^ {a + l)\ 

'' ^ Z^^ ' {d - a)\j\{a + 1 - ])r-^ ' (a + 1)! 



i=o 



^d + In v^, ,,■ ,a + In 



-(-^)^-"(:;;)D-^)^rr)(^+^) 



i=o 



Now let g(s) = s^ be a polynomial of degree /3. By HA.H , the above expression can be rewritten as 

,1/d+lx 



^a. + 1' 



(A.3) 



(1) Assume t > s. Then a + 1 > P + 1, and so by Lemma [A. II we have A°'~^^g{s) — 0. This completes the proof of the 
first part. 

(2) Let t = s - 1. By Lemma IATT] we have A°'+^g{s) = (a + 1)! and so by jPTal 

^a + l' ^ 5 ' 

This completes the proof of the second part. ■ 

Appendix B. Miscellaneous auxiliary results 

Theorem B.l (Rouche's theorem). Let the polynomial q{z) G C[z] be a sum q{z) = p{z) + e(z). Let zq be a simple zero of 
p{z). If there exists p G M+ such that 

\p{z)\ > \e{z)\ Vz e dBp{zo) 

then q{z) has a simple zero mside Bp{zo). 

Lemma B.2. Let there be given a sequence of polynomials Pk{z) : C — ^ C and a point zo G C such that 

(1) Pk(zo) = for all k e N; 

(2) Pj^(zo) > C47 for aH fc G N and some constant C47 independent of k; 
(S) For every fixed k the following inequality holds for all z G B^^-i (zq) 

|P^'(z)| <C48fc 

where Cas, i^ o. constant independent ofk. 

Let p{k) : N -> M satisfy 

< p(fc) < min (-, ^'' 



, fc C48 fc y 
Then there exists a constant C49 independent of p{k) such that for all k, the following holds for every z G dBp(k){^o)' 

\Pk{z]\ >C49p(fc) 

Proof. Let us write the truncated Taylor expansion of Pj. around zq with remainder in Lagrange form: 

Pk (zo +p(fcK«) = Pfc(zo) +P^(zo)p(fc)e'« + ^ii^^=(fc)e='« 



for some ^ G Bp{zo). Now since p[k) < i we have 



On the other hand, 



|Bi| yCi-rpik) 



pik) < ^^^ 



Cisk 

Ciskp^k) ^ Ci7P{k) 
2-2 

Therefore -L^ > IB2I and so by taking C49 '^= ^ we have \Pkiz)\ > Cic,p{k). 

Proposition B.3. Let n G N fee given. Then for \x\ < ^^ the following estimate holds: 

(1 + xy^ ^ 1 - nx + ^—^ -x'^Rz(x) where |i?5(x)|< 

^ * 2 1 a:(?i+2J 

3 

In general, for approximation of order d, we have for |ir| < "t,. 

(1 + .)-" ^l-nx+ ■:h^^.^ + ... + (-1)^ nx---x{n + d-l) ^, ^ ^_^^,+, ^12i_^^^+^,'i+i R, (,) 



where 

\R6ix]\ < 



, _ (n+d+1) 
^ d+2 ^ 



Proof. Standard majorization of the Taylor series tail by a geometric series. ■ 

Appendix C. Explicit construction op a bump 

In this appendix we present an explicit construction of the bump function which we used in our numerical simulations. We 
also derive an explicit bound for the size of its Fourier coefficients, to be used in the proof of localization accuracy. 

Let there be given two parameters t and E with 2E > t, together with the point £ G K. Our goal is to build a function 
9 ~ 5s,i(^l f ) which satisfies the following conditions: 

(Gl) g = Oforir^ [^-B^i + E]; 

(G2) 5 = lforxe [$-!,$+ |]i 

(G3) g e C°°(R); 

(G4) the Fourier coefficients of g decay as rapidly as possible. 

The idea is to take a standard C°° moUifier, scale it and convolve with a box function. 

We define two new parameters: the scaling factor s and the width of the box r. Note that our construction implies r > 2s, 
because otherwise the result of the convolution will not have a flat segment in the middle. 

Let us therefore take the standard C°° moUifier 

-1/(1-0,2) ^^^ |_^| ^ ^ 

otherwise 
and scale it between — s and 5 for some s > 0: 

where 

A = -/ *(-jdcc=/ *(y) d J/ ~ 0.443994 




Now we take a box function centered at £ , having width r 



1 for - f < X - £ < f 



otherwise 
Finally we convolve the two and get a smooth bump: 

1 r^+i fx-t\ 

g = gv,s{_x; £) = brix; £) * ms{x) = — - / * dt 

5A Jj-J V s / 

The new parameters s,r should be compatible with the original E,t. In particular, we want to have a strip of width t in the 

center, and the extent of the whole bump should not exceed E. Therefore we have the following compatibility conditions: 

t r 
s+ - < - 

2 2 (C.l) 

r 
2s+ - < E 
2 

The function g so constructed clearly satisfies conditions (Gl) — (G3) above. Let us now maximize the decay of its Fourier 
coefficients. By definition: 



27r5A 



Ck(g) = ^--r e-'^U <!'{ \dt}dx 



?- 



Notice first that ^(z) is zero outside the region —1 < z < 1, therefore we can make the change of variables z — ^t — ir,t— ^t 
and rewrite the integral as 

c,(9) = ;A r e^'^^l-) I r'^ e-'^^dt Sdz 



2-7rsA 



i-i 




-0.5 0.5 

(a) The standard mollifier ^(x) 



(b) The final bump for J = 

Figure C.l. The construction of the bump 



So now the two integrals are completely separated. Explicit calculation gives 



„ — ikt J ^ V / 



e-'^Mt 



Now we scale back: z = sy and obtain the explicit formula 









(C.2) 



27rAfc 

Finally we would like to determine the optimal values for 5 and r so that \ck {g)\ decrease as rapidly as possible with k — ^ 00. 
First note that since ^' G C°°, then for every a > 1 there exists a constant C50 (a) such that 

icfc(*)i <CBo-ifcr" 



The formula (|C.2p suggests that we should take 5 to be as large as possible. Applying the conditions (jC.ip we get that the 
following values maximize s: 

3 V 2 



r*^-(E + t) 



We have thus proved the following result: 



Theorem C.l. Given E,t with 2E > t and a point f, let 5Bt((r;£) be the bump constructed above. Then it satisfies the 
conditions (Gl) — (G4) such that for every a > 1 there exists a constant C51 = C51 (a) such that for all k £ N 



\<=k{gE,t)\ <Cbi •(2B-t)-°fc 



—a 1, — l-a 



Appendix D. Initial estimates via Prony's method 

In this appendix we present a rigorous proof that the Eckhoff 's method of order zero produces sufficiently accurate estimates 
of the jump locations {$j}, to be used in Algorithm 15. II Denote aij = e~'^3. For d = 0, the system 1)2. 2p becomes 

K 
J2 ^oj'^j « 27r (ik) CkU) (D.l) 



ID. It is a well-known system of equations which is sometimes called the Prony system ([33 ) or Sylvester-Ramanujan system 
([28]). The original method of solution (due to Baron de Prony, 32 ) is to exploit the following fact. 

Lemma D.l. The sequence {m^} satisfies the recurrence relation with constant coefficients 

K 
y mk+iqi = 



where {gi} are the coefficients of the polynomial 



Q(z) = f](z-a;,) = ^giz* 



i=i 



Proof. We have Q{uj) = for all j = 1, . . . , K. Therefore 



i—O j—1 J=l 



i=0 



i=0 



J = l 



Corollary D.2. Let qx = 1 for normalization. Then for all k ^ N the coefficient vector {gi}j_Q is the solution of the linear 
system 



(D.2) 



mk 


rrik+i 


mk+K-i 




90 




mk+K 


mk+i 


mk+2 


mk+K 


X 


91 




rrik+K+l 


mk+K-i 


mk+K ■ 


'"^k+2K-2 


y 


9K-1 




mk+2K-i 



def 

After this preparation, we can now describe the algorithm for obtaining initial estimates using Prony's method. Recall that 
our a-priori bounds are given by Ji, J2, J3 and T - see Section [5] 

Algorithm D.3. Let us be given the first M + 2K — 1 Fourier coefficients Ck{f) of a function f with K unknown discontinuities 
{fj}"'^^, which IS continuously differentiable between these discontinuities. Denote the magnitudes of the jumps by {Aj')^_.^. 

(1) Calculate the sequence 

Tk = 2-K{ik)ck{f) 



(2) Solve the system 



ru 


TM+l 


TM+K-l 


ru+i 


rM+2 


riA+K 


fM+K-l 


tm+k 


rM+2K-2 










def~ 





90 
91 



9K-1 



rM+K 
'T'M+K+l 

rM+2K-l 



(D.3) 



(3) Take the estimated < ojj \ to be the roots of the polynomial 

K-l 
Q{z) = z"^ + XI ^^-'^ 



and then set 



-argwj 



Now we would like to analyze the accuracy of Algorithm [D3] First, we need to estimate the error in solving the system (ID. 311 
We use standard result from numerical linear algebra. 

Lemma D.4. Consider the linear system Ax = b and let xo be the exact solution. Let this system be perturbed: 

{A + A^) ir = 6 + A6 

31 



and let xo + Ax denote the exact solution of this perturbed system. Denote 
\Ax\\ .. IIA^II .. IIA6II 



5x : 



^ol 



5A: 



5b: 



K = ||j4||||j4 ^11 (condition number) 



for some vector norm \\ ■ \\ and the induced matrix norm. Then 

K 



5x < 



1- K-SA 



{5A + Sb) 



(D.4) 



Proof. See e.g. 



Consider IID.311 . The error in the right-hand side is given by 115.111 . Therefore we now need to estimate the condition number of 
the matrix H^. Although all the entries are bounded, it may still happeqj that k(Hic) is unbounded. Fortunately, this is not 
the case. To see this, we are going to factorize H^ into a component which depends on fc, and a component which doesn't. 

Lemma D. 5. Let V = V ((i, . . . ,(k) denote the Vandermonde matrix on the nodes {uj}, i.e. 



1 


1 


1 


wi 


W2 


UK 


K-1 

1 


.f- . 


■ <-' 



Then for allk efi 

Hk = V X diag{>lo,ja)*} X V'^ 

Proof. Direct computation from the definitions l|D.l[l and (|D.2p . ■ 

Corollary D.6. For all k ef>i 

Remark D.7. k (V) is well-studied in e.g. [18 . It essentially depends on the minimal distance between the nodes. In particular: 

n 

1 



v'- 



n 



max 
i<i<K i- i- \ojj — Ui\ 



Lemma D.8. There exist constants Cs2,I^l2 such that for alii = 0,l,...,d and for all k > K'12— t^k(V) 

'.'1 



TVs 
J? 






Proof. In the context of Lemma [0.4 1 our original system is H^q = m l|D.2[l and the perturbed system is Jf^q = m l)D.3[l . Note 
that |mfc| > Ji ■ C53 for some C53. Rrom previous considerations we therefore have 

HH.-H.II T 1 

\\Hk\\ - Ji k 

T 1 

(5m < C55 ■ -r ■ T 

J\ k 

Jl 



We would like to estimate 5q according to HP. 411 . If 



k > 2C54 -liiV) 

^-N^-' Jl Jl 



def 



then K (Hk) SHk < 4 and so 



^ -' Jl Jl K 



def 
-C52 



Consider for instance Hj, : 



[: 1-1] 



We can finally estimate the accuracy of Algorithm ID. 31 To shorten notation, let 

Theorem D.9. For every < Q < 1 there exist constants C56 (a) ,^"13 (a, F) such that for every k > K13 Alaorithm \D.3\ 
reconstructs the locations of the jumps with accuracy 

\ii-I]\ <C56-F-k"-^ 

Proof. We have shown that the perturbation Q — Q has coefficients of magnitude k~^. We will use the same reasoning as in 
Section|4]in order to estimate the quantity lojj — i^j\ - which is the perturbation of the roots of Q (z). Let \z\ < 2. Since the 
coefficients of Q {z) do not depend on k, we can obviously find constants C57 and Css such that 

(1) \Q'{a,,)\ yCsr 

(2) \Q"{z)\ <C58 

By a reasoning similar to Lemma IB. 21 we conclude that there exist constants C59 < 1 , Ceo such that for every function 
< p{k) < C59 we have 

\Q{z)\> Ceopik) Vz e dBp^k) K) (D.5) 

Now let 

p{k) ^ F-k"-'^ 

If fc > (■^—) '^^ , then p (fc) < Csg and so (ID. 511 holds. On the other hand, by Lemma [D. 81 we have that 

|(Q-Q)(z)| <Cei -F-k-^ 

Now finally we require that k > (S^) ° , in which case 

|(Q-Q)(z)| <C6i -F-k-^ <C6oFk''-' < \Q{z)\ 

and therefore Q(z) has a simple zero inside Bpii^-\ (wj). 

Thus we have shown that Yojj — ojA < F ■ k"^^. Write Zj = Uj + P (fc) fc"^^ where \P (fc)| < F. Then by Taylor expansion of 

the logarithm we will have (recall \uj\ = 1) for large enough fc > K13 (a) 



\^:-^A 



arg Wj - arg Wj 



arg I — 

J 



< C56 (a) • F • fc° 
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